home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / langs / clisp.zoo / clisp / examples / numberth / folgasym.lsp < prev    next >
Encoding:
Lisp/Scheme  |  1991-10-25  |  9.1 KB  |  254 lines
  1. ;; Funktionen zur Bestimmung der asymptotischen Entwicklung von Folgen
  2. ;; Hartmut Maennel, Bruno Haible
  3. ;; August 1989
  4. ;; Version für Long-Floats 30.8.1989
  5.  
  6. (provide 'folgen-asymptotik)
  7.  
  8. ; Interpoliert eine Funktion an einer Wertemenge durch ein Polynom:
  9. ;   yfun(t) = sum(j>=0, c[j]*(xfun(t)-x0)^j) .
  10. ; stellen ist eine Sequence, die die Auswertungsstellen enthält.
  11. ; Wert ist ein Array der Koeffizienten #(c[0] c[1] ...)
  12. ; bei Entwicklung des Polynoms um die Stelle x0.
  13. (defun polynom-interpol (xfun yfun stellen &optional (x0 0))
  14.   (let* ((len (length stellen))
  15.          (xwerte (map 'array xfun stellen))
  16.          (ywerte (map 'array yfun stellen)))
  17.     ; Nun ist
  18.     ; xwerte[i]=xfun(stellen[i]), ywerte[i] = yfun(stellen[i]) für 0<=i<len.
  19.     ; Mit fun(x) := yfun(xfun^-1(x)) also
  20.     ; ywerte[i] = fun(xwerte[i]) für 0<=i<len.
  21.     (do ((i 0 (1+ i))) ((>= i len))
  22.       (do ((j (1+ i) (1+ j))) ((>= j len))
  23.         (setf (aref ywerte j)
  24.               (/ (- (aref ywerte j) (aref ywerte i))
  25.                  (- (aref xwerte j) (aref xwerte i))
  26.     ) ) )     )
  27.     ; Nun ist fun(x) = 
  28.     ; = ywerte[0]+(x-xwerte[0])*(ywerte[1]+(x-xwerte[1])*(ywerte[2]+...))
  29.     ; für alle x=xwerte[i], 0<=i<len.
  30.     (do ((i (1- len) (1- i))) ((< i 0))
  31.       ; Multipliziere das Polynom
  32.       ; ywerte[i+1]*(x-x0)^0+...+ywerte[len-1]*(x-x0)^(len-i-2) mit
  33.       ; x-xwerte[i] = (x-x0) + (x0-xwerte[i]) und addiere ywerte[i]:
  34.       (let ((c (- x0 (aref xwerte i))))
  35.         (do ((j i (1+ j)))
  36.             ((>= j (1- len)))
  37.           (incf (aref ywerte j) (* c (aref ywerte (1+ j))))
  38.     ) ) )
  39.     ; Nun ist fun(x) = ywerte[0]*(x-x0)^0+ywerte[1]*(x-x0)^1+...
  40.     ; für alle x=xwerte[i], 0<=i<len.
  41.     ywerte
  42. ) )
  43.  
  44. (setf (long-float-digits) 665) ; entspricht 200 Dezimalstellen
  45.  
  46. ;(defparameter pi (* 4 (atan 1L0)))
  47. (defparameter e (exp 1L0))
  48.  
  49. (defun myrationalize (a) ; liefert eine rationale Zahl in der Nähe von a.
  50.                          ; Kriterium: Kettenbruch-Element >10000
  51.   (when (floatp a) (setq a (rational a)))
  52.   (if (< a 0)
  53.     (- (myrationalize (- a)))
  54.     (let ((p0 1) (q0 0) (p1 (floor a)) (q1 1) next)
  55.       (setq a (- a p1))
  56.       (loop
  57.         (when (zerop a) (return))
  58.         (multiple-value-setq (next a) (floor (/ a)))
  59.         (when (> next 10000) (return))
  60.         (psetq p0 p1 p1 (+ p0 (* p1 next))
  61.                q0 q1 q1 (+ q0 (* q1 next))
  62.       ) )
  63.       (/ p1 q1)
  64. ) ) )
  65.  
  66. ; mehr Stellen ergeben größere Rechenzeit, aber auch größere Genauigkeit:
  67. ; Der relative Fehler in limes ist etwa (1/n)^(Stellenzahl-1) = 1/n^k.
  68. (defun default-k (&optional (max 100))
  69.   (cond ((>= max 400) 5)
  70.         ((>= max 100) 8)
  71.         ((>= max 20) 10)
  72.         (t (floor max 2))
  73. ) )
  74. (defun default-stellen (&optional (max 100) (k (default-k max)))
  75.   ; Liste von max/(1+i/k), i=0..k:
  76.   (do ((i 0 (1+ i))
  77.        (L nil))
  78.       (nil)
  79.     (push (round max (+ 1 (/ i k))) L)
  80.     (when (>= i k) (return L))
  81. ) )
  82.  
  83. (defun interpol-ntel (fun &optional (max 100) (k (default-k max))
  84.                                     (stellen (default-stellen max k)))
  85.   (polynom-interpol #'/ fun stellen)
  86. )
  87.  
  88. (defun limes (fun &optional (max 100) (k (default-k max)))
  89.   (aref
  90.     (interpol-ntel
  91.       #'(lambda (n) (coerce (funcall fun n) 'long-float))
  92.       max
  93.       k
  94.     )
  95.     0
  96. ) )
  97.  
  98. ; Bestimmt die Glieder c0, c1, ..., cm einer asymptotischen Entwicklung
  99. ; fun(n) = c0 + c1*n^-1 + c2*n^-2 + ... + cm*n^-m
  100. (defun asymptotik (fun m &optional (max 100) (k (default-k max)))
  101.   (let ((c (make-array (1+ m)))
  102.         (Lc (make-array (1+ m)))
  103.         (j 0)
  104.         (hauptnenner 1)) ; Hauptnenner von c[0]..c[j-1]
  105.     (loop
  106.       (setf (aref c j)
  107.         (/ (myrationalize
  108.              (* hauptnenner ; Hauptnenner erleichtert myrationalize
  109.                 (limes
  110.                   #'(lambda (n)
  111.                       (let ((y (funcall fun n)))
  112.                         (dotimes (i j) (setq y (* n (- y (aref Lc i)))))
  113.                         y
  114.                     ) )
  115.                   max
  116.                   k
  117.            ) )  )
  118.            hauptnenner
  119.       ) )
  120.       (print (aref c j))
  121.       (setq hauptnenner (lcm hauptnenner (denominator (aref c j))))
  122.       (setf (aref Lc j) (coerce (aref c j) 'long-float))
  123.       (incf j)
  124.       (when (> j m) (return))
  125.     )
  126.     c
  127. ) )
  128.  
  129.  
  130. ; Aufwand bei  max = ca. 100,  Stellenzahl  k+1 = ca. 6:
  131. ; Funktion polynom-interpol löst ein LGS mit einer Matrix, die eine
  132. ; Vandermonde-Matrix zu (2/max, ..., 1/max) ist und daher die Determinante
  133. ; prod (0<=i<j<=k, ((1+j/k)/max - (1+i/k)/max) )
  134. ; = 1!2!...k! / (k max)^(k(k+1)/2) hat. (Davon der Logarithmus ->)
  135. ; Dies verbraucht etwa  k^2/2 * (lg e + lg max) Stellen.
  136. ; Erfahrungsgemäß hat ein so gewonnener limes eine relative Genauigkeit
  137. ; von 1/max^k, also  k * lg max  genaue Stellen.
  138. ; Will man  m (= ca. 5)  Koeffizienten der asymptotischen Entwicklung aus-
  139. ; rechnen, so muß man die Werte um  m * lg max  genauer haben:
  140. ; Bei (...((Wert - c0) * max - c1) * max ...) gehen  m * lg max  Stellen
  141. ; verloren.
  142. ; Somit muß man etwa  m * lg max + k^2/2 * lg (e*max)  Stellen vorsehen und
  143. ; bekommt dann alle Limites mit  k * lg max genauen Stellen, und der
  144. ; Rechenaufwand ist:
  145. ;   m *                  m mal Limes ziehen
  146. ;   k^2 *                k^2 Rechenoperationen
  147. ;   (k^2 lg max)^2       auf Zahlen mit  k^2 lg max  Stellen
  148. ;  = m * (lg max)^2 * k^6.
  149. ; Also:
  150. ; Sofern die Folge leicht berechenbar ist (z.B. durch polynomial lineare
  151. ; Rekursion), max groß und k klein wählen (z.B. max=400 und k=5 -> 13 Stellen
  152. ; beim Limes, brauche 85 Stellen und bei m=10 zusätzlich 26 Stellen).
  153. ; Sofern nur endlich viele Folgenglieder vorliegen, k größer wählen
  154. ; (was die Rechenzeit steil in die Höhe treibt): z.B. max=50 und k=9 ->
  155. ; 15 Stellen beim Limes, brauche 96 Stellen und bei m=5 zusätzlich 8 Stellen.
  156.  
  157. ; Man beachte, daß man einen Koeffizienten c_i EXAKT haben muß, um c_i+1
  158. ; ausrechnen zu können. Das Ganze ist nämlich extrem empfindlich gegenüber
  159. ; Fehlern. Um c_i exakt (aus den ersten 14 Stellen des Dezimalbruchs) zu
  160. ; bestimmen, kann myrationalize dienen (falls c_i rational ist), oder
  161. ; der LLL-Algorithmus (falls c_i in einem algebraischen Erweiterungskörper
  162. ; der rationalen Zahlen liegt).
  163.  
  164. (require 'intmisc) ; defun-N0
  165.  
  166. #|
  167. ; Beispiel:
  168. (require 'kartenf "KARTEN4")
  169.  
  170. ; Das Originalproblem: (kartenf n)
  171.  
  172. ; Davon die Fakultätsfaktoren (4n)!/n!^4 entfernen:
  173. (defun fa (n) (/ (kartenf n) (/ (! (* 4 n)) (expt (! n) 4))))
  174.  
  175. ; Probiere (limes #'(lambda (n) (/ (fa n) (fa (1- n)))) 100),
  176. ; sollte etwa 25/64 liefern.
  177.  
  178. ; Davon die Potenzen (5/8)^2n entfernen:
  179. (defun fb (n) (/ (fa n) (expt 5/8 (* 2 n))))
  180.  
  181. ; Probiere (limes #'fb 100), sollte etwa (sqrt 125/216) liefern.
  182.  
  183. ; Dann in Long-Float umwandeln:
  184. (defun fc (n) (coerce (fb n) 'long-float))
  185.  
  186. ; Dann konstanten Faktor (5/6)^(3/2) entfernen:
  187. (defvar c00 (sqrt (coerce 125/216 'long-float)))
  188. (defun-N0 f0 (n) (/ (fc n) c00)) ; bis hierher stets nur einmal berechnen
  189.  
  190. ; Dann die asymptotische Entwicklung bestimmen:
  191. ; (asymptotik #'f4 2 100)
  192. ; liefert (1 -35/432 -815/41472) als Anfang der asymptotischen Entwicklung.
  193. |#
  194.  
  195. ; oder schrittweise:
  196. ; c0 = (myrationalize (limes #'f0)) = 1
  197. (defun f1 (n) (* n (- (f0 n) c0)))
  198. ; c1 = (myrationalize (limes #'f1))
  199. (defun f2 (n) (* n (- (f1 n) c1)))
  200. ; c2 = (myrationalize (limes #'f2))
  201. (defun f3 (n) (* n (- (f2 n) c2)))
  202. ; c3 = (myrationalize (limes #'f3))
  203. (defun f4 (n) (* n (- (f3 n) c3)))
  204. ; c4 = (myrationalize (limes #'f4))
  205. (defun f5 (n) (* n (- (f4 n) c4)))
  206. ; c5 = (myrationalize (limes #'f5))
  207. (defun f6 (n) (* n (- (f5 n) c5)))
  208. ; c6 = (myrationalize (limes #'f6))
  209. (defun f7 (n) (* n (- (f6 n) c6)))
  210. ; c7 = (myrationalize (limes #'f7))
  211. (defun f8 (n) (* n (- (f7 n) c7)))
  212. ; c8 = (myrationalize (limes #'f8))
  213. (defun f9 (n) (* n (- (f8 n) c8)))
  214. ; c9 = (myrationalize (limes #'f9))
  215. (defun f10 (n) (* n (- (f9 n) c9)))
  216. ; c10 = (myrationalize (limes #'f10))
  217. (defun f11 (n) (* n (- (f10 n) c10)))
  218. ; c11 = (myrationalize (limes #'f11))
  219. (defun f12 (n) (* n (- (f11 n) c11)))
  220. ; c12 = (myrationalize (limes #'f12))
  221. (defun f13 (n) (* n (- (f12 n) c12)))
  222. ; c13 = (myrationalize (limes #'f13))
  223. (defun f14 (n) (* n (- (f13 n) c13)))
  224. ; c14 = (myrationalize (limes #'f14))
  225. (defun f15 (n) (* n (- (f14 n) c14)))
  226. ; c15 = (myrationalize (limes #'f15))
  227. (defun f16 (n) (* n (- (f15 n) c15)))
  228. ; c16 = (myrationalize (limes #'f16))
  229. (defun f17 (n) (* n (- (f16 n) c16)))
  230. ; c17 = (myrationalize (limes #'f17))
  231. (defun f18 (n) (* n (- (f17 n) c17)))
  232.  
  233. #|; zusammengefaßt:
  234. (defun f1 (n) (* n (- (f0 n) c0)))
  235. (defun f2 (n) (* n (- (f1 n) c1)))
  236. (defun f3 (n) (* n (- (f2 n) c2)))
  237. (defun f4 (n) (* n (- (f3 n) c3)))
  238. (defun f5 (n) (* n (- (f4 n) c4)))
  239. (defun f6 (n) (* n (- (f5 n) c5)))
  240. (defun f7 (n) (* n (- (f6 n) c6)))
  241. (defun f8 (n) (* n (- (f7 n) c7)))
  242. (defun f9 (n) (* n (- (f8 n) c8)))
  243. (defun f10 (n) (* n (- (f9 n) c9)))
  244. (defun f11 (n) (* n (- (f10 n) c10)))
  245. (defun f12 (n) (* n (- (f11 n) c11)))
  246. (defun f13 (n) (* n (- (f12 n) c12)))
  247. (defun f14 (n) (* n (- (f13 n) c13)))
  248. (defun f15 (n) (* n (- (f14 n) c14)))
  249. (defun f16 (n) (* n (- (f15 n) c15)))
  250. (defun f17 (n) (* n (- (f16 n) c16)))
  251. (defun f18 (n) (* n (- (f17 n) c17)))
  252. |#
  253.  
  254.